Method and system for numerical simulation of fluid flow

ABSTRACT

A method of numerical simulation of fluid flow past a body. The physics of stationary flow is formulated in terms of a potential functional of at most four scalar fields: three Clebsch scalar fields and a density field. The physics of non-stationary flow is formulated in terms of an action that includes an action integral of a Lagrangian functional of the same four scalar fields and an initial and a final integral of a function of the same four scalar fields. The values of the fields are varied to extremize the potential functional or the action, under the constraint of appropriate boundary conditions. Potential functionals and Lagrangian functionals for compressible and incompressible flows with zero velocity normal to the surface of the body, and for potential flows, are described.

FIELD AND BACKGROUND OF THE INVENTION

[0001] The present invention relates to the numerical simulation of fluid flow and, more particularly, to a numerical method of simulating barotropic fluid flow past a body by minimizing a potential functional.

[0002] Barotropic fluid flow is described by the Euler equations: $\begin{matrix} {{\frac{\partial\overset{\rightarrow}{v}}{\partial t} + {\overset{\rightarrow}{v} \cdot {\overset{\rightarrow}{\nabla}\overset{\rightarrow}{v}}} + {\overset{\rightarrow}{\nabla}\left( {h + \Phi} \right)}} = 0} & (1) \end{matrix}$

[0003] and of the continuity equation: $\begin{matrix} {{\frac{\partial\rho}{\partial t} + {\overset{\rightarrow}{\nabla}{\cdot \left( {\rho \quad \overset{\rightarrow}{v}} \right)}}} = 0} & (2) \end{matrix}$

[0004] where {right arrow over (ν)}(x^(k),t) is the velocity vector field, Φ(x^(k),t ) is a potential relating to an external force, ρ(x^(k),t) is the scalar density field, and h, the specific enthalpy, is a function of ρ through the equation of state of the fluid. {right arrow over (ν)}, Φ and ρare functions of position ^(k)=(x,y,z) in space and of time t. In the special case of a stationary flow with small viscosity, the time derivatives vanish:

{right arrow over (ν)}·{right arrow over (∇)}{right arrow over (ν)}={right arrow over (∇)}(h+Φ)=0;{right arrow over (∇)}×(ρ{right arrow over (ν)})=0  (3)

[0005] Only in trivial cases can these equations be solved analytically. In cases of practical interest, for example, in aerodynamics and hydrodynamics, these equations must be solved numerically. The most common way to solve these equations is by integrating them. The possible numerical instabilities associated with numerical integration are well known and need not be detailed here. A numerical solution based on a variational principle would be inherently numerically stable. Ideally, such a numerical solution would involve finding the minima of a potential functional with respect to the associated variables. In the case of fluid flow, the variables are four scalar fields: ρ and the three Cartesian components of {right arrow over (ν)}. Heretofore, no variational formulation of barotropic fluid flow has represented the solutions of equations (1) and (2) or of equations (3) as the values of {right arrow over (ν)} and ρ, or of any other set of only four scalar fields, that minimize a potential functional. The smallest number of scalar fields obtained heretofore was seven, by R. L. Seliger and G. B. Witham (Proc. Roy. Soc. London, Vol. A305, p. 1, 1968). Because equations (1)-(3) depend on only four scalar fields, it should be possible to obtain a variational formulation of barotropic fluid flow in terms of only four scalar fields.

[0006] There is thus a widely recognized need for, and it would be highly advantageous to have, a method of finding numerical solutions of the equations of fluid flow by minimizing a potential functional of four scalar fields.

SUMMARY OF THE INVENTION

[0007] According to the present invention there is provided a numerical method of simulating fluid flow past a body having a boundary, including the steps of: (a) formulating the fluid flow in terms of a potential functional of at most four firmament scalar fields and at most two scalar variables; and (b) extremizing the potential functional with respect to the at most four fundamental scalar fields and with respect to the at most two scalar variables.

[0008] According to the present invention there is provided a system for numerical simulation of fluid flow, including: (a) a software module including a plurality of instructions for computing, from discrete representations of at most four fundamental scalar fields, a potential functional representative of the fluid flow, and for varying the discrete representations to extremize the potential functional; (b) a processor for executing the instructions; and (c) a memory for storing the discrete representations.

[0009] According to the present invention there is provided a numerical method of simulating fluid flow, including the steps of: (a) formulating the fluid flow in terms of a potential functional of at most four fundamental scalar fields and at most two scalar variables; and (b) extremizing the potential functional with respect to the at most four fundamental scalar fields and with respect to the at most two scalar variables.

[0010] According to the present invention there is provided a numerical method of simulating fluid flow past a body having a boundary, including the steps of: (a) formulating the fluid flow in terms of an action which is the sum of: (i) an action integral, from an initial time t₀ to a final time t₁, of a Lagrangian functional L of at most four fundamental scalar fields, (ii) a spatial integral of a function of the at most four fundamental scalar fields at the final time t₁, and (iii) a negative of a spatial integral of the function of the at most four fundamental scalar fields at the initial time t₀; and (b) extremizing the action with respect to the at most four fundamental scalar fields.

[0011] According to the present invention there is provided a numerical method of simulating fluid flow including the steps of: (a) formulating the fluid flow in terms of an action which is the sum of: (i) an action integral, from an initial time t₀ to a final time t₁, of a Lagrangian functional L of at most four fundamental scalar fields, (ii) a spatial integral of a function of the at most four fundamental scalar fields at the final time t₁, and (iii) a negative of a spatial integral of the function of the at most four fundamental scalar fields at the initial time to; and (b) extremizing the action with respect to the at most four fundamental scalar fields.

[0012] According to the present invention there is provided a system for numerical simulation of fluid flow, including: (a) a software module including a plurality of instructions for computing, from discrete representations of at most four fundamental scalar fields, an action representative of the fluid flow, and for varying the discrete representations to extremize the action; (b) a processor for executing the instructions; and (c) a memory for storing the discrete representations.

[0013] It is shown in the Appendix that stationary barotropic fluid flow can be formulated as the extremization of a potential fimctional of at most four scalar fields: the Clebsch scalar fields α, β and ν and the density field ρ; and also of at most two scalar variables β₁ and ν₁l. Specifically, the functional is minimized with respect to α, β and ν and maximized with respect to β₁ and ν₁. The scalar fields are functions of the position vector x^(k). The velocity vector field {right arrow over (ν)} is related to the Clebsch scalar fields through {right arrow over (ν)}=α{right arrow over (∇)}β+{right arrow over (∇)}ν.

[0014] The potential functionals are as follows:

[0015] If the component of {right arrow over (ν)}normal to the boundary of the region of fluid flow is zero, or if the density p on the boundary is zero, the potential functional is $\begin{matrix} {{\int_{V}{\left\lbrack {{\frac{1}{2}\quad \left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} + {ɛ(\rho)} + \Phi} \right\rbrack \rho {^{3}x}}} - {\beta_{1}\left( {{\int_{V}{{\alpha\rho}\quad {^{3}x}}} - {\overset{—}{\alpha}\quad M_{0}}} \right)} - {v_{1}\left( {{\int_{V}{\rho \quad {^{3}x}}} - M_{0}} \right)}} & (4) \end{matrix}$

[0016] (equations 36 and 42 of the Appendix) for compressible flows and $\begin{matrix} {{\int_{V}{\left\lbrack {{\frac{1}{2}\left( {{\alpha \quad {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} + \Phi} \right\rbrack \rho_{0}{^{3}x}}} - {\beta_{1}\left( {{\int_{V}{{\alpha\rho}_{0}\quad {^{3}x}}} - {\overset{\_}{\alpha}\quad M_{0}}} \right)}} & (5) \end{matrix}$

[0017] (equations 45 and 46 of the Appendix) for incompressible flows. ε(ρ)is the specific internal energy of the fluid, in the case of a compressible flow, and is related to the specific enthalpy by $\begin{matrix} {h = \frac{\partial\left( {\rho \quad ɛ} \right)}{\partial\rho}} & (6) \end{matrix}$

[0018] {overscore (α)} and M₀ are constants. The integrals are over the volume V occupied by the fluid. ρ₀ is the constant density, in the case of an incompressible flow. The functional is minimized with respect to α, β and ν, and maximized with respect to β₁. In the case of a compressible flow, the functional also is minimized with respect to ρ and maximized with respect to ν₁l.

[0019] In the case of potential flows, for which {right arrow over (∇)}×{right arrow over (ν)}=0, the potential functional is $\begin{matrix} {{\int_{V}{\left\lbrack {{\frac{1}{2}\left( {\overset{\rightarrow}{\nabla}v} \right)^{2}} + {ɛ(\rho)} + \Phi} \right\rbrack \rho \quad {^{3}x}}} - {v_{1}\left( {{\int_{V}{\rho \quad {^{3}x}}} - M_{0}} \right)}} & (7) \end{matrix}$

[0020] (equations 82 and 83 of the Appendix) for compressible flows and $\begin{matrix} {\int_{V}{\left\lbrack {{\frac{1}{2}\left( {\overset{\rightarrow}{\nabla}v} \right)^{2}} + \Phi} \right\rbrack \rho_{0}\quad {^{3}x}}} & (8) \end{matrix}$

[0021] (equation 86 of the Appendix) for incompressible flows. ε(92 ) is the specific internal energy of the fluid, in the case of a compressible flow. The integrals are over the volume occupied by the fluid. M₀ is a constant. ρ₀ is the constant density, in the case of an incompressible flow. The functional is minimized with respect to ν. In the case of a compressible flow, the functional also is minimized with respect to ρ and maximized with respect to ν₁.

[0022] It also is shown in the Appendix (see equation 93 therein) that non-stationary barotropic fluid flow can be formulated as the maximization of an action which is the sum of three integrals: an action integral of a Lagrangian functional L of three Clebsch scalar fields (α, βand ν and the density field ρfrom an initial time t₀ to a final time t₁: $\begin{matrix} {\int_{t_{0}}^{t_{1}}{L\quad {t}}} & (9) \end{matrix}$

[0023] an initial integral over the four scalar fields at time t₀:

−∫(α(x ^(k) ,t ₀)β(x ^(k) ,t ₀)+ν(x ^(k) ,t ₀))ρ(x ^(k) ,t ₀)d ³ x  (10)

[0024] and a similar final integral over the four scalar fields at time t1:

(∫α(x ^(k) ,t ₁)β(x ^(k) ,t ₁)+ν(x ^(k) ,t ₁))ρ(x ^(k,) t, ₁)d ³ x  (11)

[0025] The spatial integrals are over the volume V occupied by the fluid. Note that the integrand of the initial integral is the negative of the integrand of the final integral.

[0026] If the component of {right arrow over (ν)} normal to the boundary of the region of fluid flow is zero, or if the density p on the boundary is zero, the Lagrangian functional is: $\begin{matrix} {\int_{V}{\left\lbrack {{- \left( {{\alpha \frac{\partial\beta}{\partial t}} + \frac{\partial v}{\partial t}} \right)} - {\frac{1}{2}\left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} - {ɛ(\rho)} - \Phi} \right\rbrack \rho \quad {^{3}x}}} & (12) \end{matrix}$

[0027] (equation 11 of the Appendix) for compressible flow and $\begin{matrix} {\int_{V}{\left\lbrack {{- \left( {{\alpha \frac{\partial\beta}{\partial t}} + \frac{\partial v}{\partial t}} \right)} - {\frac{1}{2}\left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} - \Phi} \right\rbrack \rho_{0}\quad {^{3}x}}} & (13) \end{matrix}$

[0028] (equation 15 of the Appendix) for incompressible flows. As in the case of stationary flows, ε(ρ) is the specific internal energy of the fluid, in the case of a compressible flow; the integrals are over the volume V occupied by the fluid; and ρ₀ is the constant density, in the case of an incompressible flow. The maximization is with respect to α, β and ν, and also with respect to ρ in the case of a compressible flow.

[0029] The Clebsch scalar fields (α, β and ν and the density ρ are referred to herein as “fundamental” scalar fields because only these fields are varied explicitly to minimize the potential functionals. The other scalar fields that appear in the potential functionals are functions of these fundamental scalar fields. For example, ε is a function of ρ, via equation (6).

[0030] In addition to being inherently stable, the present invention, as applied to stationary flows, is faster than the prior art methods and is easier to use for complex geometries.

A NOTE ON NOTATION

[0031] In the Appendix, the symbols α, β and ν are used to represent Clebsch variables that are functions of both space x^(k) and time t. In the case of stationary flows, the space and time dependences of the Clebsch variables can be separated as follows:

α={circumflex over (α)}(x ^(k))  (14)

β={circumflex over (β)}(x ^(k))−β₁ t  (15)

ν={circumflex over (ν)}(x ^(k))−ν₁ t  (16)

[0032] In potential functionals for stationary flows as derived in the Appendix, the Clebsch scalar fields are these functions {circumflex over (α)}, {circumflex over (β)} and {circumflex over (ν)} of x^(k) only; and in the course of the derivations, the dependence on time drops out, leaving β₁ and ν₁ as scalar variables to be varied, along with the values of the scalar fields {circumflex over (α)}, {circumflex over (β)} and {circumflex over (ν)}, to extremize the potential functionals. For simplicity elsewhere herein, when the potential functionals of stationary flows are discussed, the Clebsch scalar fields that enter into these potential functionals are denoted by α, β and ν, without carets.

[0033] The symbol that is used herein for velocity, a lowercase italic “V” (ν), and the symbol that is used herein for the third Clebsch scalar field, a lowercase Greek “nu” (ν), are typographically similar. To avoid confusion, the reader should bear in mind that the symbol for velocity always appears with an arrow above it ({right arrow over (ν)}) when referring to the vectorial velocity field, or with the subscripts x or y when referring to Cartesian components of the velocity field; and that the symbol for the third Clebsch scalar field appears either unadorned (ν) when referring to the third Clebsch scalar field itself, or with the subscript 1 when referring to a constant that is related to the third Clebsch scalar field.

BRIEF DESCRIPTION OF THE DRAWINGS

[0034] The invention is herein described, by way of example only, with reference to the accompanying drawings, wherein:

[0035]FIG. 1A is a contour plot of the x-component of the velocity of flow past a cylinder, computed analytically;

[0036]FIG. 1B is a contour plot of the x-component of the velocity of flow past a cylinder, computed numerically according to the present invention;

[0037]FIG. 2A is a contour plot of the y-component of the velocity of flow past a cylinder, computed analytically;

[0038]FIG. 2B is a contour plot of the y-component of the velocity of flow past a cylinder, computed numerically according to the present invention;

[0039]FIG. 3A is a contour plot of the magnitude of the velocity of flow past a cylinder, computed analytically;

[0040]FIG. 3B is a contour plot of the magnitude of the velocity of flow past a cylinder, computed numerically according to the present invention;

[0041]FIGS. 4A and 4B are plots of sections of the Clebsch scalar field ν, computed both analytically and numerically according to the present invention;

[0042]FIG. 5 is a high level block diagram of a system of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0043] The present invention is of a method of computing the velocity and pressure of a barotropic fluid flowing past a body.

[0044] The principles and operation of computational fluid dynamics according to the present invention may be better understood with reference to the drawings and the accompanying description.

[0045] For numerical computation, the continuous scalar fields α, β, ν and ρ, and functions thereof such as ε(ρ) are represented by their values at discrete points on regular grids. In the Appendix, two interleaving rectangular grids are described for 2-dimensional computations, one grid for βand ν,and the other grid for α, ρ, Φ and h. It will be clear to those skilled in the art how to extend these grids to three dimensions. The grids occupy the portion of space occupied by the fluid. In other words, values of the scalar fields are not specified inside the body.

[0046] Initial values of the Clebsch fields α, β and ν and, for compressible flow, of the density ρ, are specified at the grid points, and these values are varied to minimize the appropriate potential functional. The actual fields of interest are the velocity {circumflex over (ν)} and the density ρ. It is shown in the Appendix how to obtain the initial values of the Clebsch fields from the initial values of the velocity. The final values of the velocity are obtained from the final values of the Clebsch fields from {right arrow over (ν)}=α{right arrow over (∇)}β+{right arrow over (∇)}ν({right arrow over (ν)}=ν in the case of potential flows).

[0047] The grids on which the scalar fields are represented are finite in size. For numerical simulation of flow past a finite body, the grids are made large enough so that the values of the velocity at the outer boundaries of the grids are the same as they would be in the absence of the body; and the velocity is fixed at those values.

[0048] In the case of the component of {right arrow over (ν)} normal to the boundary of the region of fluid flow being zero, it is not necessary to specify in advance the values of the tangential components of the velocity at the boundary of the body. Simply finding the values of α, β, ν and β₁ (and also ρ and ν₁ for compressible flows), at grid points outside the body, that minimize the potential functional (4) or (5) automatically produces a velocity field {right arrow over (ν)}=α{right arrow over (∇)}β+{right arrow over (∇)}ν that has a vanishing component normal to the body.

[0049] Similarly, in the case of potential flows, merely varying the values of the scalar fields only on grid points outside the body, to minimize the appropriate potential functional, yields the desired solution.

[0050] In the case of the density ρ on the boundary being zero, the boundary condition that is specified in advance at the grid points closest to the boundary of the body is precisely that: ρ=0.

[0051] If the component of {right arrow over (ν)} normal to the boundary of the region of fluid flow is zero, or if the density ρ on the boundary is zero, the constants {overscore (α)} and M₀ also must be specified. M₀ is the volume integral of the initial estimate of ρ. {overscore (α)} is the volume integral of the initial estimate of αρ, divided by M₀.

[0052] Referring now to the drawings, FIGS. 1-4 compare the analytic solution and the numerical solution of the present invention for incompressible fluid flow past a circular cylinder. For an incompressible fluid flowing past a cylinder of radius α, if the flow is in the x-direction and the velocity of the fluid is asymptotically U at x=±∞, the analytic solution is: $\begin{matrix} {v = {{Ux}\left( {1 + \frac{a^{2}}{x^{2} + y^{2}}} \right)}} & (17) \\ {v_{x} = {U\left( {1 - {a^{2}\left( \frac{x^{2} - y^{2}}{\left( {x^{2} + y^{2}} \right)^{2}} \right)}} \right)}} & (18) \\ {v_{y} = {{- 2}U\frac{a^{2}}{\left( {x^{2} + y^{2}} \right)^{2}}{xy}}} & (19) \\ {{\overset{\rightarrow}{v}} = \sqrt{v_{x}^{2} + v_{y}^{2}}} & (20) \end{matrix}$

[0053] The numerical solution was obtained by minimizing potential (8) with respect to the Clebsch scalar field ν, represented by its values on a rectangular grid of 14641 points (121×121) with a spacing of 0.1 meters. The cylinder had a radius of 0.6 meters and the center of the cylinder was coincident with the center point of the grid. The asymptotic value U of the velocity was 1 m/sec. In accordance with equation (75) of the Appendix, the initial value of v was a linear scalar field whose gradient was U in the x-direction and zero in the y-direction. Φ was taken to be zero. Note that the minimum of potential (8) is independent of ρ₀, so that the numerical problem was one of varying v to minimize the absolute value of its gradient, subject to the boundary conditions on the surface of the cylinder.

[0054]FIG. 1A is a contour plot of the analytical values of ν_(x). FIG. 1B is a contour plot of the corresponding numerical values of ν_(x). FIG. 2A is a contour plot of the analytical values of ν_(y). FIG. 2B is a contour plot of the corresponding numerical values of ν_(y). FIG. 3A is a contour plot of the analytic values of velocity magnitude |{right arrow over (ν)}|. FIG. 3B is a contour plot of the corresponding numerical values of |{right arrow over (ν)}|. The axis labels are grid point indices.

[0055]FIG. 4A is a plot of ν along the line y=0 m from x=−3 m to x=3 m. FIG. 4B is a plot of ν along the line x=0.8 m from y=−3 m to y=3 m. Note that the points −0.6 m ≦×≦0.6 m in FIG. 4A are inside the cylinder. In both FIGS. 4A and 4B, the solid line represents analytic values of ν and the discrete points represent numerical values of ν.

[0056] The numerical computations for non-stationary flows are similar to the numerical computations for stationary flows, with the addition of one more dimension (time) to the grids. The time that a system requires to settle down from an initial non-stationary condition to a final stationary condition usually can be estimated, by dimensional analysis, as the ratio of a characteristic length to a characteristic velocity. This time is taken as t₁ -t₀. The values of the Clebsch scalar fields and of the density field are specified at time t₀ to define the problem. The same values define the starting values of the scalar fields at all the other grid points, and those values are varied to minimize the action. Alternatively, the stationary flow technique is used to determine the final condition of the system, and the starting values of the scalar fields at intermediate times are determined by linear interpolation. Note that in general it is necessary to specify and fix the values of each of four scalar fields (three scalar fields for incompressible flows) at one particular time (t₀, t₁, or any time in between) to define the time-dependent problem. For example, the time-dependent problem can be defined by specifying initial and final values of β and ν throughout space. The time-dependent problem defined in this way is solved by varying stand ρ at all times, and β and ν at intermediate times, to minimize the action.

[0057]FIG. 5 is a high level block diagram of a system 10 for numerical simulation of fluid flow according to the present invention. System 10 includes a processor 12, a random access memory 14 and a set of input/output devices, such as a keyboard, a floppy disk drive, a printer and a video monitor, represented by I/O block 16. Memory 14 includes an instruction storage area 18 and a data storage area 20. Within instruction storage area 18 is a software module 22 including a set of instructions which, when executed by processor 12, enable processor 12 to simulate fluid flow by the method of the present invention.

[0058] Using the appropriate input device 16 (typically a floppy disk drive), source code of software module 22, in a suitable high level language, for numerical simulation of fluid flow according to the present invention, i.e., by extremizing a potential with respect to at most four scalar fields, or by maximizing an action with respect to at most four scalar fields, is loaded into instruction storage area 18. Selecting a suitable language for the instructions of software module 22 is easily done by one ordinarily skilled in the art. The language selected should be compatible with the hardware of system 10, including processor 12, and with the operating system of system 10. Examples of suitable languages include but are not limited to compiled languages such as FORTRAN, C and C++. If a compiled language is selected, a suitable compiler is loaded into instruction storage area 18. Following the instructions of the compiler, processor 12 turns the source code into machine-language instructions, which also are stored in instruction storage area 18 and which also constitute a portion of software module 22. Using the appropriate input device 16 (typically a keyboard), the parameters of the fluid flow problem are entered, and are stored in data storage area 20. Following the machine-language instructions, processor 12 stores the initial values of discrete representations of up to four scalar fields in data storage area 20 and varies those values to extremize the potential, selected from among potentials (4), (5), (6), (7) and (8), that is appropriate to a stationary fluid flow problem, or to maximize the action defined by integrals (9) through (11) and the appropriate Lagrangian functional, selected from among Lagrangian functionals (12) and (13), that is appropriate to a non-stationary fluid flow problem. The results of the extremization are displayed at video monitor 16 or printed on printer 16, preferably in graphical form as in FIGS. 1-4.

[0059] It will be appreciated that any of the relevant numerical integration methods known in the art for, for example finite element methods and multi-grid methods, may be used to evaluate the potential functionals of the present invention. Furthermore, as is well known in the art, the coordinate system used need not be Cartesian, but may be curvilinear.

[0060] While the invention has been described with respect to a limited number of embodiments, it will be appreciated that many variations, modifications and other applications of the invention may be made. 

What is claimed is:
 1. A numerical method of simulating fluid flow past a body having a boundary, comprising the steps of: (a) formulating the fluid flow in terms of a potential functional of at most four fundamental scalar fields and at most two scalar variables; and (b) extremizing said potential functional with respect to said at most four fundamental scalar fields and with respect to said at most two scalar variables.
 2. The method of claim 1, wherein said at most four fundamental scalar fields include three Clebsch scalar fields α, β and ν.
 3. The method of claim 1, wherein said at most four fundamental scalar fields include a density field.
 4. The method of claim 1, wherein said at most four fundamental scalar fields include: (a) three Clebsch scalar fields α, β and ν and (b) a density field ρ.
 5. The method of claim 4, wherein said potential functional is ${{\int_{V}{\left\lbrack {{\frac{1}{2}\left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} + {ɛ(\rho)} + \Phi} \right\rbrack \rho \quad d^{3}x}} - {\beta_{1}\left( {{\int_{V}{{\alpha\rho}\quad d^{3}x}} - {\overset{\_}{\alpha}M_{0}}} \right)} - {v_{1}\left( {{\int_{V}{\rho \quad d^{3}x}} - M_{0}} \right)}};$

β₁ and ν₁ being said two scalar variables; ε being a specific internal energy, Φ being a potential of an external force, α{overscore (α)} and M₀ being constants and V being a volume of integration exclusive of the body; said extremizing being effected with respect to α, β, ν, ρ, β₁ and ν₁.
 6. The method of claim 1, wherein the fluid flow is formulated in terms of a potential functional of three Clebsch scalar fields α, βand ν; and wherein said extremizing is effected with respect to α, β and ν and with respect to only one said scalar variable.
 7. The method of claim 6, wherein said potential functional is ${\int_{V}{\left\lbrack {{\frac{1}{2}\left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} + \Phi} \right\rbrack \rho_{0}d^{3}x}} - {\beta_{1}\left( {{\int_{V}{{\alpha\rho}_{0}d^{3}x}} - {\overset{\_}{\alpha}M_{0}}} \right)}$

β₁ being said only one scalar variable, Φ being a potential of an external force, ρ₀, {overscore (α)} and M₀ being constants and V being a volume of integration exclusive of the body.
 8. The method of claim 1, wherein said fluid flow is formulated in terms of a potential functional of a single Clebsch scalar field ν and a density ρ; said extremizing being effected only with respect to ν and ρ.
 9. The method of claim 8, wherein said potential functional is ${{\int_{V}{\left\lbrack {{\frac{1}{2}\left( {\overset{\rightarrow}{\nabla}v} \right)^{2}} + {ɛ(\rho)} + \Phi} \right\rbrack \rho \quad d^{3}x}} - {v_{1}\left( {{\int_{V}{\rho \quad d^{3}x}} - M_{0}} \right)}};$

wherein ε is a specific internal energy, Φ is a potential of an external force, M₀ is a constant and V is a volume of integration exclusive of the body.
 10. The method of claim 1, wherein said fluid flow is formulated in terms of a potential functional of a single Clebsch variable ν; said extremizing being effected only with respect to ν.
 11. The method of claim 10, wherein said potential functional is $\int_{V}{\left\lbrack {{\frac{1}{2}\left( {\overset{\rightarrow}{\nabla}v} \right)^{2}} + \Phi} \right\rbrack \rho_{0}d^{3}x}$

wherein Φ is a potential of an external force, ρ₀ is a constant and V is a volume of integration exclusive of the body.
 12. A system for numerical simulation of fluid flow, comprising: (a) a software module including a plurality of instructions for computing, from discrete representations of at most four fundamental scalar fields, a potential functional representative of the fluid flow, and for varying said discrete representations to extremize said potential functional; (b) a processor for executing said instructions; and (c) a memory for storing said discrete representations.
 13. A numerical method of simulating fluid flow, comprising the steps of: (a) formulating the fluid flow in terms of a potential functional of at most four fundamental scalar fields and at most two scalar variables; and (b) extremizing said potential functional with respect to said at most four fundamental scalar fields and with respect to said at most two scalar variables.
 14. A numerical method of simulating fluid flow past a body having a boundary, comprising the steps of. (a) formulating the fluid flow in terms of an action which is the sum of: (i) an action integral, from an initial time t₀ to a final time t₁, of a Lagrangian functional L of at most four fundamental scalar fields, (ii) a spatial integral of a function of said at most four fundamental scalar fields at said final time t₁, and (iii) a negative of a spatial integral of said function of said at most four fundamental scalar fields at said initial time t₀; and (b) extremizing said action with respect to said at most four fundamental scalar fields.
 15. The method of claim 14, wherein said at most four fundamental scalar fields include three Clebsch scalar fields α, β and ν.
 16. The method of claim 14, wherein said at most four fundamental scalar fields include a density field.
 17. The method of claim 14, wherein said at most four fundamental scalar fields include: (a) three Clebsch scalar fields α, β and ν and (b) a density field ρ.
 18. The method of claim 17, wherein said Lagrangian functional is $\int_{V}{\left\lbrack {{- \left( {{\alpha \frac{\partial\beta}{\partial t}} + \frac{\partial v}{\partial t}} \right)} - {\frac{1}{2}\left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} - {ɛ(\rho)} - \Phi} \right\rbrack \rho \quad d^{3}x}$

wherein ε is a specific internal energy, Φ is a potential of an external force, and V is a volume of integration exclusive of the body.
 19. The method of claim 14, wherein the fluid flow is formulated in terms of a Lagrangian functional of three Clebsch scalar fields α, β and ν; and wherein said extremizing of said action is effected with respect to said α, β and ν.
 20. The method of claim 19, wherein said Lagrangian functional is $\int_{V}{\left\lbrack {{- \left( {{\alpha \frac{\partial\beta}{\partial t}} + \frac{\partial v}{\partial t}} \right)} - {\frac{1}{2}\left( {{\alpha {\overset{\rightarrow}{\nabla}\beta}} + {\overset{\rightarrow}{\nabla}v}} \right)^{2}} - \Phi} \right\rbrack \rho \quad d^{3}x}$

wherein Φ is a potential of an external force, ρ₀ is a constant and V is a volume of integration exclusive of the body.
 21. A numerical method of simulating fluid flow comprising the steps of: (a) formulating the fluid flow in terms of an action which is the sum of: (i) an action integral, from an initial time t₀ to a final time t₁, of a Lagrangian functional L of at most four fundamental scalar fields, (ii) a spatial integral of a function of said at most four fundamental scalar fields at said final time t₁, and (iii) a negative of a spatial integral of said function of said at most four fundamental scalar fields at said initial time t₀; and (b) extremizing said action with respect to said at most four fundamental scalar fields.
 22. A system for numerical simulation of fluid flow, comprising: (a) a software module including a plurality of instructions for computing, from discrete representations of at most four fundamental scalar fields, an action representative of the fluid flow, and for varying said discrete representations to extremize said action; (b) a processor for executing said instructions; and (c) a memory for storing said discrete representations. 